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ABSTRACT 

We investigate the means by which cold gas can accrete onto Milky Way mass galaxies 
from a hot corona of gas, using a new smoothed particle hydrodynamics code, 'SPHS'. 
We find that the 'cold clumps' seen in many classic SPH simulations in the literature 
are not present in our SPHS simulations. Instead, cold gas condenses from the halo 
along filaments that form at the intersection of supernovae-driven bubbles from previ- 
ous phases of star formation. This positive feedback feeds cold gas to the galactic disc 
directly, fuelling further star formation. The resulting galaxies in the SPH and SPHS 
simulations differ greatly in their morphology, gas phase diagrams, and stellar con- 
tent. We show that the classic SPH cold clumps owe to a numerical thermal instability 
caused by an inability for cold gas to mix in the hot halo. The improved treatment of 
mixing in SPHS suppresses this instability leading to a dramatically different physical 
outcome. In our highest resolution SPHS simulation, we find that the cold filaments 
break up into bound clumps that form stars. The filaments are overdense by a factor 
of 10-100 compared to the surrounding gas, suggesting that the fragmentation results 
from a physical non-linear instability driven by the overdensity. This 'fragmenting 
filament' mode of disc growth has important implications for galaxy formation, in 
particular the role of star formation in bringing cold gas into disc galaxies. 
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1 INTRODUCTION 

The star formation rate (SFR) of the Universe has 
been rapidly falling since a redshift of ~ 2 — 3 (e. g., 
iLillv et al.. 199dMada"u" et al. . 199dHippelein et al.. 2003t ). 
dominated by the most massive galaxies and galaxy 
groups. By contrast, the Milky Way - similarly to 
other disc galaxies - has been forming stars at a near- 
continuous rate for the past ^ 8 Gyr sin ce z ^ 1 (e.g., 
iNoh fc Scalo. 199CiRocha-Pinto et al.. 2000l ). The observed 
SFR of ~ 1 — 3 Mq yr~^ is hard to explain given the amount 
of cold gas present in the disc today; the Milky Way must 
have continuously accreted cold gas at a rate of ^ 1 M q 
yr~^ over this time (|Fraternali fc Tomassetti. 20121 ). 
Merger- driven acc retion appears to account for just 
~ 0.1 Mq yr"^ jSancisi et al.. 20081 ). and there is to 
date no evi dence of a low reds hift 'cold flow' accretion 
mode (e.g. IStewart et al.. 201ll ). This has led to two 
main solutions in the literature. The first is recycling of 
gas fro m existing stars in the disc through stellar wind s 
(e.g., [Roberts. 196jSandage. 198dKennicutt et al.. 1993 ). 



iLeitner fc Kravtsov (201ll ) estimate that this recycled gas 
could contribute at least half of the global SFR for a galaxy 
of Milky Way mass at low redshift {z ^ 0.5). The second 
is accretion from a massive hot halo, or corona, of gas 
surrounding the Galaxy. Such a hot halo has not yet been 
directly observed, but several indirect li nes of evidence exis t: 
observations of X-ray emitting gas (|Gupta et al.. 20121 ): 
absorption along sight lines to quasars (e.g. , 
IWiUiams et al.. 200dFang et al.. 200dKacprzak et al.. 2008 ): 

pulsar dispersion measures (e-g-, 

[Anderson fc Bregman. 201ClGaensler et al.. 20081 ): a 
significant Galactic baryon deficiency when com- 
pared to the universal baryon fra ction (e.g., 
iFukugita fc Peebles. 200dNicastro et al.. 20081 ) ; and ev- 
idence of ram pressure st ripping of the Magellanic stream 
jMastropietro et al.. 2005|) and other nearby dwarf galax- 
ies ( Grcevich fc Putman. 20091 ). These studies give a 
hot halo mass of > 5 X lO'^ - 1.5 X 10^° Mq assummg a 
iNavarro et al. ( 19961 ) (NFW) profile or > 4 x 10^° assumin g 



a flattened power-law profile ([Anderson fc Bregman. 201ll ) 
Similar results are seen in other disc galaxies (e 



iRead fc Trentham. 200dMo et al.. 200dSancisi et al.. 200aAnderson fc Bregr 
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While it is likely that hot gaseous coronae sur- 
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round disc galaxies, it is not clear how the gas cools 
and condenses out of these coronae and onto the disc 
to form stars. One particular mode of cold gas supply 
that has been seen in a number of numerical simulations 



briefly review the SPHS algorithm and 'classic' SPH. We de- 
scribe our initial conditions, and present our implementation 
of radiative cooling, star formation and stellar feedback, and 

our treatment of a cent ral supermassive black hole (SMBH). 

200 1iI€taffilac)aiHl8l(&l31w(2liM)g i^tmani^alts 2<i^ we discuss in 
Section [S] Finally, in Section [SI we present our conclusions. 



2 METHODS 
2.1 SPHS 

The SPHS me thod is described in det ail in 
iRead et al. f2010bl ) and iRead fc Havfield (20121 ). We 
give a brief summary of the main equations here. The 
discretised hydrodynamic eq uations of motion are chosen 
to minimise force errors as in lRead et al. (2010bl ): 



Pi ^^mjWij{\rij\,hi) 



ISommer-Larsen. 200dKaufmann et al.. 200dKaufmann et al. 
but was initially proposed bv lNulsen (IQSsi l. is that of direct 
cooling from the halo via thermal instability. However , 
do ubt has been cast on this picture bv lMalagoli et al." (19871 ) 
& iBinnev et al. (20091 ) who show that hot haloes are lin- 
early sta ble to density pertu rbations, malting direct cooling 
unlikely. Ijoung et al. (2013 ) extend this treatment to the 
non-linear regime with dedicated numerical simulations, 
finding that while a runaway process of cooling and collapse 
is possible, it can only occur for overdensities of ;i 10 — 20 
with respect to the local background density. These re- 
sults suggest that if gas is to cool from the hot coronae 
then some mechanism is requir ed to excite a thermal 
instability. iMarinacci et al. (201ll ) suggest a mechanism 
whereby a galactic fountain seeds metal-rich gas into the 
metal-poor hot haloes, giving rise to a thermal instability 
that causes cold gas to rain do wn onto the disc in the form 
of ~ 10^ Mq clouds (see also iFraternali fc Binnev. 20o3 ). 
This model provides an excellent fit to both the kine- 
matics and sp atial distribution of warm HI gas in the 
Milky Way (|Marasco et al.. 20121 ). although it can- 
not account for the high- velocity clouds (HVCs) (e. g., 
ISembach et al. 200jrripp et al.. 200jCollins et al.. 2005h . 
Alternative models include cooling stripped gas from dwarf 
galajcies or warm clouds at the disc-corona interface ( e.g., 
iPutman et al.. 200dHeitsch fc Putman. 200dPeek. 20091 ). 

The 'direct cooling' mode mentioned above 
is seen regularly in 'classicQ smoothed par- 

ticle hydrodynamics (SPH) simulations (e.g., 

ISommer-Larsen. 200dKaufmann et al. 200dKaufmann et al.. 200l^4.\iMmsmr&.iM)1^mi^n^leml^h(Myi ). 

lim W(\r-r'\,h) = 5(|r-r'|) 



~dt 



, Pip 3 



P. = A^p] 



(1) 



(2) 



(3) 



where rrii is the mass of particle i; Wij = 
^[Wij{hi) + Wij{hj)]; and W is a symmetric kernel 
that obeys the normalisation condition: 



W{\r-r'\,h)d\' 



where a thermal instability leads to a sudden and widespread 
condensation of gas from the halo in the form of cold, dense 
clumps. If correct, no special mechanism would be re- 
quired to explain the continued star formation of disc 
galaxies over the past ~ 8 Gyrs. However, 'classic' SPH 
is known to exhibit an artificial surface tension that 
inhibits mixing of different gaseous phases, leading to poor 
performance on a variety of hydrodynamic test problems 
lAgertz et al.. 2007|Price. 200dWadslev et al. 200dRead et al.. 
In recent work, some of the present authors have developed 
a new 'flavour' of SPH - SPHS - that solves these problems, 
giving excellent performance and numerical convergenc e 
on a wide range of test problems (|Read fc Havfield. 20ld ). 
In this paper, we use SPHS to revisit the problem of 
thermal instabilities in hot gaseous coronae. Our primary 
goals are to determine whether the instabilities seen in 
the classic SPH simulations are physical or numerical; 
and under what circumstances cold gas can condense 
out of a hot halo to fuel star formation in disc galaxies. 
The former goal is further motivated by the absence of 
such cold clumps in the hot haloes of galaxy formation 
simul ations using adaptiv e mesh refinement (AMR) codes 
(e.g., lAgertz et al.. 20091) or the recent mo ving Voronoi 
mesh code AREPO (|Vogelsberger et al.. 201ll ). 

This paper is organised as follows. In Section [2] we 



^ for a careful definition of 'classic' SPH see Section 12.21 



(4) 



(5) 



and Tij — r j — is the vector position of the particle relative 
to the centre of the kernel. 

We use a variable s moothing length hi as in 
ISpringel fc Hernguist (2002al ) that is adjusted to obey the 
following constraint equation: 



2010al ). 



Air 



hiUi — Nn ; with m — Wij (6) 



where Nn is the typical neighbour number (the number of 
particles inside the smoothing kernel, W). The above con- 
straint equation gives fixed mass inside the kernel if particle 
masses are all equal. For the full SPHS runs (see Section 
m we use the 'HOCT4' kernel with 442 neighbours as this 
gives significantly improved force accuracy and convergence 



gives signmcantly improved torce accuracy ; 
(IRead et al. . 2010ttlead fc Havfield. 20121) : 



w — — < 

h^ ^ 
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< a:: < K 
K, < X < j3 
P < X <a 
X <1 

otherwise 



with TV = 6.515, P 
and K = 0.214. 



-2.15, Q = 0.981, a = 0.75, ^ 
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In addition to the above equations of motion, numer- 
ical dissipation is switched on if particles are converging. 
This avoids multivalued fluid quantities occurring at the 
point of convergent flow. Without such dissipation, the re- 
sulting multivalued pressures drive waves through the fluid 
that propagate large numerical errors and spoil convergence. 
The switch is given by: 



fe?|V(V-Vi)l 



V ■ Vi < 
otherwise 



(8) 

where aioc,i describes the amount of dissipation for a given 
particle in the range [0, Ofmax = 1]; Cj is the sound speed 
of particle i; and ris = 0.05 is a 'noise' parameter that de- 
termines the magnitude of velocity fluctuations that trigger 
the switch. Equation [8] turns on dissipation if V • Vi < 
(convergent flow) and if the magnitude of the spatial deriva- 
tive of V • Vi is large as compared to the local divergence 
(i.e., if the flow is going to converge). The key advantage as 
compared to most other switches in the literature is that it 
acts as an early warning system, switching on before large 
numerical errors propag ate throughout the fluid (see also 
ICuUen fc Dehnen. 20ld V The second derivatives of the ve- 
locity field are calculated usin g high order polynomia l gra- 
dient estimators descri bed in iMaron fc Howes (20031 ) and 
iRead fc Havfield (20121 '). We use the above switch to turn 
on dissipation in all advected fluid quantities - i.e., the mo- 
mentum (artificial viscosity) and entropy (artificial thermal 
conductivity). Once the trajectories are no longer converg- 
ing, the the dissipation parameter decays back to zero on a 
timescale ~ hi/ci. The dissipation equations are fully con - 
servative and described in detail in lRead fc Havfield (2012 ). 
The only free parameter is amax that describes the rate of 
dissipation that occurs when particle trajectories attempt 
to cross. Note that by construction, the numerical dissipa- 
tion occurs at the resolution limit of the simulation and so 
convergence is independent of «max. 

In addition to the hydrodynamical modifications, the 
SPHS method also includes an improved timestepping 
algorithm for strong shocks. This was adapted from 
ISaitoh fc Makino (20091 ). who found that the evolution of a 
Sedov- Taylor blast wave (or similar) is captured incorrectly 
when the gas in the shock is on a very different timestep to 
the gas it is impinging upoijfl. In extreme cases, the parti- 
cles sitting ahead of the shock may not 'feel' the particles in 
the blast wave, as the stationary particles are on such long 
timeste ps the shock has passed by before their next force 
update. ISaitoh fc Makino (20091 ) recommend restricting the 
timesteps between neighbours to be at most a factor of 4 in 
order to alleviate this problem; this is what we do also in 
SPHS. 

The SPHS algorithm was incorpo rated into th e 
N-body/hydrodynamical code GADGET-3 (|Springel. 2005h . 
which was used for all the simulations in this paper. 



^ Wc would like to thank Frazer Pearce & Stewart Muldrew for 
providing the first version of the timestepping code patch that wo 
use. 



2.2 Classic SPH 

Throughout this paper, we will present comparisons 
with 'classic' SPH. We define this to be the fully 
conservative 'entropy' form of SPH described in 



ISpringel fc Hernguist (2002lf ) . However, many pr oduc- 
tion SP H codes - e.g. Gasoline (|Wadslev et al.. 20041 ). and 
Hydra (jCouchman et al.. 19951 ) - are sufRciently similar 
to warrant also being described as 'classic' SPH. The 
discretised Euler equations are the same as in SPHS, but 
with the momentum equation replaced by: 



Pi Pj 



(9) 



where the function fi is a correction factor that ensures en- 
ergy conservation for varying smoothing lengths: 



.A = 1 + 



hi dpi 
3pi dhi 



(10) 



We do not use the above conservative momentum equation 
in SPHS since it leads to larger force errors with only a 
modest improvement in energy conservation (at least when 
ap plied to galaxy and gala xy cluster formation sim ulations; 
see lRead et al. (2010bl ) and lPead fc Havfield (20ll ) for fur- 
ther details). 

Unlike in SPHS, there is no dissipation switching and 
Q = Qmax = const. = 1 always. There is also no dissipation 
in entropy; the only numerical dissipation applied is the ar- 
tificial viscosity. This prevents multivalued momenta from 
occurring, but not multivalued entropy or pressure (e.g., 
iRead et al.. 2010bh . 

For most simulations presented in this paper, we use 
a standard cubic spline (CS) kernel with 96 neighbours 
- both for SPH and SPHS. This is somewhat larger 
than usually employed but is necessary for the high or- 
d er gradients required for the dissipation switch in SPHS 
((Read fc Havfield. 20121 ). By using the same neighbour num- 
ber for both hydrodynamic techniques, we ensure like spa- 
tial resolution. However, at this neighb our number, the CS 
kernel is prone to a pairing instability (|Read et al.. 2010al ) 
that could in principle seed numerical instabilities in a hot 
corona. We explicitly check that this is not the case by run- 
ning one of our classic SPH simulations with a more stan- 
dard 42 neighbours, and with the HOCT4 kernel with 442 
neighbours (that is manifestly stable to pairing). In all cases, 
we recover the result seen in the literature that the hot coro- 
nae breaks up into many cold clumps. 



2.3 Initial conditions 

O ur setup is similar to the ' cooling gaseous halo' m odel 
of iKaufmann et al. (20061 ) and iKaufmann et al. (20071 ) . For 
our initial condition, we employ a live dark matter (DM) 
halo of coUisionless particles together with a gaseous 
halo of SPH particles, with each component relaxed for 
many dynamical times with an adiabatic equation of 
state (EQS) in order to remove Poisson noise. The relax- 
ation step is performed separately for each resolution and 
kernel-neighbour number combination used. Both the DM 
and gas distributions follow a Dehnen-McLaughlin model 
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l|Dehnen fc McLaughlin. 20051 ) of the form 

1 

p(r) cx 



(r/rj/^) [l + (r/rf ^)] 



(11) 



with the gas profile normalised to contain a mass of 1.5 x 
lO" M0 and the DM to contain a mass of 1.5 x 10^^ M©. 
The scale radius for the haloes is = 40 kpc, and both are 
truncated at a virial radius rt = 200 kpc. 

The gas is initially (after relaxation) in hydrostatic equi- 
librium, with a temperature profile given by the relation 



T{r) 



fimp 



kB Pgas(r) 



GM(r) 



dr 



(12) 



where ^ is the mean molecular weight, fca is the Boltz- 
mann constant, pgas(r) is the radial gas density profile 
and M(r) is the enclosed mass of both components within 
a radius r. The gas is given a velocity field whereby 
the specific angular momentum profile foll ows a power-law 
l|Bullock et al.. 200]|Kaufmann et al.. 20071 ) such that 



(13) 



and normalised by 
by 



rotation parameter A — 0.038, defined 



A = 



.7 gas I -E^dn 



|l/2 



GM. 



3/2 



(14) 



where -Edm and Mdm are the total energy and mass of the 
DM halo. This normalisation implicitly assumes negligible 
angular momentum transport between the DM halo an d the 
gcis, an assumption justified in lKaufmann et al. (20071 ). The 
rotational velocity field is implemented about the z-axis. 

Relevant parameters for the simulations are given in Ta- 
ble [1] The DM haloes were constructed using variable par- 
ticle masses in order to mi nimise computatio nal expense, as 
described in Section 3.1 of lCole et al. (20111 ). Gravitational 
softening lengths for the coUisionless N-body components 
(DM & stars) were fixed, with the values given in Table 
[B (tdm, which also applies to the softening lengths of any 
stars formed during the simulation). For the gas, we em- 
ployed adaptive gravitational softening lengths which were 
set equal to the gas smoothing length at all times so as to 
ensure there was no artificial bias towards pressure support 
or gravitational collap se in Jeans-unstable regions, as per 
iBate fc Burkert (19971 ). We did not employ the conservative 
correction terms as suggested in [Pricc fc Monaghan (20071) : 
we will explore these in future work. To explicitly check 
that our variable softening does not affect our results, we 
re-ran the SPH-96-resl and SPHS-96-resl simulations with 
fixed gravitational softening in the gas, set to the minimum 
smoothing length in the latter two simulations. Our results 
were largely unchanged in terms of the SPH/SPHS compar- 
ison. 



2.4 Radiative cooling 

For all of the simulations, we employ radiative cool- 
ing in the gas with a cooling floor of Tfloor ~ 
100 K. We use a toy cooling curve that follows 
iKatz et al. (19961) above lO" K. assu ming primordial abun- 
dance, and Mashchenko et al. (20081 ) below 10* K assuming 



Solar abundance. This crudely models a low metallicity cool- 
ing halo that rapidly reaches ~ Solar abundance in cooling 
star forming regions. We will consider a more realistic cool- 
ing curve that self-consistently responds to the injection of 
metals from star forming regions in future work. 

In addition to a cooling floor, we also prevent gas from 
cooling to the point at which the Jeans mass for gravita- 
tional collapse becomes unresolved. A given SPH/SPHS 
simulation has a mass resolution equal to Afros ~ A'^rosWigas, 
where nigas is the mass of a gas particle and A'rcs is the 
number of gas particles that c onstitutes a resovable m ass, 
i.e., a single resolution element. iBate fc Burkert" (19971 ) find 
this to be approximately 2A''noigh, where A^ncigh is the num- 
ber of neighbours of a gas particle; this number is somewhat 
of a rule of thumb and depends on the resolving scale of 
the particular kernel employed. For all of the tests we use 
Arcs = 128. The resolving scales of our CS and H0CT4 ker- 
nels (see Section r2.1|l are such that this number is reasonable 
(for more detail on the resolving power of these kernels see 
iRead et al.. 2010£tlead fc Havfield. 201 jPehnen fc Alv. 20121 ). 

Our 'dynamic' cooling fioor effectively ensures that the 
Jeans mass is always resolved within our simulation. For 
a given mass element Mrcs, we can write a Jeans density, 
namely. 



PJ = 



/ TvkT 



(15) 



which manifests in the simulation as a polytropic equa- 
tion of state P = yl(s)p*''^. Gas is not allowed to 
collapse (for a given temperature) to densities higher 
than given by equation 1151 and we identify gas that lies 
on the polytrope as 'star- forming', allowing it to form 
stars above a fixed density threshold (see Section 12. 5|) . 
We emphasise that the polytrope is not physical, but 
a purely numerical device to prevent star-forming gas 
from becoming unresolved in our simulations (see, e.g., 
iTruelove et al.. 1997lvan de Voort et al.. 201llDubois et al.. 20121 ). 



2.5 Star formation and feedback 

Star formation (SF) in our simulations is modelled in a 'sub- 
grid' fashion according to observational constraints on the 
SF rate and efficiency. We allow gas that lies on the poly- 
trope to form stars (N-body particles) above a fixed density 
threshold, with an e fficiency of 0.1 as per observations (e.g., 
iLada fc Lada. 20031 ) of giant molecular clouds (GMCs). The 
star formation rate follows the Sch midt volume density law 
for star formation (|Schmidt. 1959h . namely 



PSPR a: Pgas 



(16) 



which is implemented in the simulation by employing the 
dynamical time as the relevant SF timescale, i.e. 



dp. 
At 



= V 



Pga 
idyi 



(17) 



where rj is the SF efficiency. 

In all the simulations, we include feedback from super- 
novae (SNe), which takes the form of an injection of thermal 
energy from the star particle into nearby gas particles. Due 
to our finite resolution, each star particle is not an individ- 
ual star but is instead representative of a stellar distribution. 
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Table 1. Parameters for each simulation. Two resolutions were run for each numerical method. A'^ is the total number of particles for 
each species, rrigas is the (constant) mass of each gas particle and rridm gives the range of masses for the dark matter particles. A/rcs is 
the 'resolvable' gas mass, defined in Section 12.41 h refers to the smoothing length of the gas while e refers to the gravitational softening 
length. For all simulations, e = h for the gas particles. 



ID 


iVgas 


JVdm 


TVlgas (M0) 


md 


m (Mq) 


Aires (Mq) 


/imin (kpc) 


hmax (kpc) 


Edm (kpc) 


SPH-96-resO 


1.5 X 10''^ 


1.8 X 10^ 


1 X lO'' 


1 - 


14 X 10'' 


1.28 X 10* 


0.1 


1790 


0.45 


SPHS-96-resO 


1.5 X 10'"^ 


1.8 X 10^ 


1 X 10'' 


1 - 


14 X 106 


1.28 X 10** 


0.08 


295 


0.45 


SPH-96-resl 


7.5 X 10^ 


9 X 10^ 


2 X 10^ 


2 - 


29 X IQS 


2.56 X 10^ 


0.03 


310 


0.2 


SPHS-96-resl 


7.5 X 10^ 


9 X 10^ 


2 X 10^ 


2 - 


29 X 10^ 


2.56 X 10^ 


0.03 


230 


0.2 


SPHS-442-res2 


3.75 X 10'' 


4.5 X 10^ 


4 X 10-* 


4- 


59 X Iffi 


5.12 X lO'' 


0.008 


200 


0.03 



and so we integrate over a Salpeter IMF (|Salpeter. 19551 ) be- 
tween 8 Mq and 100 Mq in order to determine the number 
of SNe that should be feeding back at any particular time. 
We treat only type II SNe, with the energy of each super- 
nova event set to Bsn = 10^^ ergs. To determine the time 
of injection we use the standard relation: 



(18) 



tMS ^ f M_ 

to ~ VMo, 

to determine an approximate main sequence time Ims after 
the initial formation, at which the energy injection is imple- 
mented as a delta function in time. We couple the feedback 
energy to a given mass of gas, rather than just the neigh- 
bours of the star particle - this mass is set to be our resolv- 
able mass Mros (see Section[23]). In doing so, we ensure that 
both the star formation and the feedback is resolved. The 
actual injection of energy is convolved with a CS kernel to 
ensure a smoothing of the thermal coupling over the mass 
scale receiving the supernova feedback. 

Th polytropic pressure floor, together with the fixed SF 
density threshold, is crucial to our sub-grid modelling strat- 
egy for star formation; it means that the gas that is able 
to form stars is (i) undergoing Jeans collapse and (ii) at 
the typical density of GMC star-forming regions. By con- 
straining 'star-forming gas' to lie on the polytrope we are 
ensuring that star formation occurs in collapsing gas re- 
gions, rather than just in regions of high density. We note 
that there has been some recent evidence that GMCs and 
GMC complexes may not need to be globally gravitationally 
bound for regions within them to undergo star formation 
(jPobbs et al.. 2011f ) but that even in this picture the dense 
sub-structures actually forming stars are both (i) in a state 
of gravitational collapse and (ii) at a density of ^ 10^ atoms 
cm"^. 

A further aspect of our polytropic pressure floor is to 
maintain consistency with the parameters of the SF sub-grid 
model. These are tailored to star-forming gas at approxi- 
mately a GMC mass scale, and employ the corresponding 
global observational constraints, rather than including any 
of the flner details of the star formation process (e.g., metal- 
line cooling, non-equilibrium chemistry, etc.). To allow gas 
to collapse to higher and higher densities as it cools would 
lead to a situation whereby not only would this collapse be- 
come dominated by numerical effects but also reach a state 
that was inconsistent with the global constraints of our sub- 
grid prescription. We plan to explore our sub-grid modelling 
strategy in this context in a forthcoming paper (Hobbs et 



al., in preparation), with particular importance placed on 
convergence with increasing resolution. In the current paper 
we do not expect convergence of the star formation proper- 
ties owing to the fact that the polytrope line shifts to lower 
mass scales as the resolution is increased; this however is 
not a problem for the comparison we make between the two 
numerical methods, since this is done at identical resolution. 
We note that such a 'non-convergent' approach is similar to 
what is typically used in other galaxy formation sub-grid 
models in the literature. 

2.6 SMBH sub-grid model 

We model the SMBH at the cent re of the halo as a sink par- 
ticle (see, e.g.. lBate et al.. 19951 ). The accretion of gas onto 
the SMBH is implemented using an 'accretion radius' ap- 
proach; gas that makes it inside a distance of race =0.1 kpc 
from the black hole is removed from the simulation, while 
its mass is then added to the SMBH mass Mbh- We start the 
SMBH as a seed of negligible mass since even the accretion 
of a single gas particle grows it to a mass of ~ 10'* — 10® Mq 
depending on the resolution (see Table 1). As a result of 
this mass being negligible before any accretion occurs, it was 
necessary to tie the black hole to the centre of the poten- 
tial in order to stop it wandering excessive distances during 
the early phases of accretion. We therefore re-position the 
SMBH on the centre of the potential at every timestep. 

For the simulations presented in this paper we do not 
model feedback from the SMBH, using it as a sink for ac- 
creting gas only. Our prescription for SMBH growth in the 
present paper is somewhat crude, but we stress that this 
is because our focus is on the differences between the SPH 
and SPHS methods in terms of thermal instabilities in the 
hot halo gas; keeping our SMBH sub-grid model simple al- 
lows us to present a cleaner numerical test between the two 
hydrodynamic methods. We plan to investigate the role of 
black hole feeding and feedback in similar galaxy formation 
simulations with a more advanced sub-grid model in future 
papers. 



3 RESULTS 

We run four main simulations, two in classic SPH, and two 
in SPHS each at two different resolutions (see Table [!}. The 
global properties between the different runs are similar due 
to the identical initial conditions; as the gaseous halo cools 



6 Alexander Hobbs, Justin Read, Chris Power, David Cole 



it loses hydrostatic support and begins to collapse, with the 
small net rotation velocity about the z-axis causing a disc 
to form initially in the x — y plane. We start by describing 
the overall behaviour of the simulations, and then discuss the 
differences in the results between the two numerical methods 
and any resolution-dependant departures from this overall 
behaviour. Our fiducial comparison is made between the two 
higher resolution runs, SPH-96-resl and SPHS-96-resl. 

As the gaseous halo collapses, it falls towards the centre 
of the computational domain. Within the first ~ 10 million 
years, some of the gas that has reached the central ~ 1 
kpc undergoes a starburst, while the rest is accreted by the 
SMBH. The feedback associated with this star formation 
event is strong enough to drive a near-spherical overpres- 
surised bubble outwards into the infalling gas, temporarily 
evacuating the central region of most of its gas and halting 
the infall within the expanding cavity that is created. Star 
formation largely ceases at this point, while the gas that 
has formed into a small-scale disc in the central kiloparsec 
continues to accrete at a significantly lower rate. The ex- 
panding shell sweeps up mass and slows down, eventually 
allowing further infall to the central regions, and in most 
cases, subsequent starbursts and feedback events. 

While the initial starburst was largely spherical, driv- 
ing a very symmetric bubble into the surrounding gas, the 
subsequent feedback events occur with a greater degree of 
asymmetry, and so hot bubbles are generated at a variety of 
different orientations at different times. In these later events 
gas is able to flow in more readily along the intersections be- 
tween bubbles, where it is able to build up the disc in the 
central few kiloparsecs. 

3.1 Differences between SPH and SPHS 

3.1.1 Lower resolution (resO) 

We start by comparing the low resolution runs, SPH-96-resO 
and SPHS-96-resO. These are shown at a time oi t = 2 Gyr 
in Figure [T] The initial evolution is similar, with some of 
the gas that reaches the central kiloparsec being accreted 
by the SMBH, and some undergoing a starburst that drives 
a hot bubble into the gaseous halo. In the SPHS case, this 
initial feedback is slightly stronger, with the gas that has 
received the thermal 'kick' reaching higher temperatures - 
a maximum of 5 x 10* K vs. 8 x 10^ K - and expanding 
slightly faster than in the SPH case. Immediately after the 
feedback event, the accretion rate drops by several orders of 
magnitude (see Figure [S]). The expanding shell slows down 
as it sweeps up mass from the gaseous halo, and it is at 
this point that the two simulations differ; in the SPH run 
the gas is able to make its way back inside the central few 
kiloparsecs, where it cools and forms stars, giving rise to a 
second starburst at ~ 1.2 Gyr. In the SPHS case, however, 
the gas has not collapsed back into the central regions by 
the end of the simulation at t = 2 Gyr, and so there is no 
associated subsequent starburst event by this time. 

We can understand some of these differences by looking 
at the star formation history between the two low resolution 
runs, which can be seen in Figure [5] along with the accretion 
history. The SMBH accretion rate largely follows the star 
formation rate, and so in the SPH run there are 2-3 major 
starburst and accretion events that take place to later times 



after the initial one. In the SPHS case the initial starburst 
lasts for longer - rather than undergoing a sudden event and 
then dropping to zero, as it does in SPH, the star formation 
rate shows a more sustained period of activity for the first 
~ 0.5 Gyr. This is the reason for the lack of the second star- 
burst in the SPHS run - the sustained initial star formation 
and feedback has driven the gas to higher temperatures and 
lower densities and so it takes a much longer time to cool 
and condense back into the central regions, failing to do so 
by the end of the simulation at t — 2 Gyr. 

Crucially, even at this low resolution we can see that by 
far the main difference between the SPH and SPHS simu- 
lations is the presence of overdense clumps that have con- 
densed from the ambient halo gas between ~ 10 — 50 kpc. 
These clumps are not present in the SPHS run, and we go 
into more detail on this particular result for the higher res- 
olution case in Section fS.l.SI 



3.1.2 Higher resolution (resl) 

The higher resolution runs constitute our fiducial compari- 
son between the two methods. Gas surface density projec- 
tions at t = 1.2 Gyrs are shown in Figure^ Once again the 
initial starburst and feedback events are similar, although 
the initial feedback is slightly stronger in SPHS than in 
SPH, with the kicked gas reaching a maximum of 10** K 
vs. 2 X 10^ K and again expanding faster. At this resolu- 
tion the peak of the star formation rate for this starburst is 
slightly higher in SPHS. Due to the stronger feedback, more 
of the gas is evacuated in the initial overpressurised bubble, 
and takes longer to make its way back inside the expand- 
ing shell. The accretion rate in SPHS is therefore reduced 
for the period immediately following the first starburst (be- 
tween ~ 0.2 — 0.7 Gyr) compared to SPH (see FigurelS)). 

Also at this higher resolution there are subsequent ma- 
jor starbursts in both the SPH and SPHS runs. The second 
large starburst occurs at ~ 1 Gyr, and takes place slightly 
earlier in SPH than in SPHS. This time the second SPH star- 
burst is more powerful than the SPHS one. We attribute this 
to the fact that immediately before the second starburst the 
SPH run saw multiple condensing gas clumps entering the 
central few kiloparsecs, which was not seen in SPHS. The 
expanding bubbles in the second starburst also differ notice- 
ably between the two methods - in SPH the boundaries of 
the cavities are thinner whereas in SPHS these outer 'walls' 
are more mixed in with the surrounding gas (see Figure [2]). 

The star formation and accretion histories for the higher 
resolution runs can be seen in Figure (6] Most noticeable is 
the presence of three peaks in the SFR to later times in the 
SPHS run, compared to two peaks in SPH. Once again the 
accretion rate largely follows the SFR, although not pre- 
cisely, and in particular we notice that whenever there is a 
strong dip in the SFR, the corresponding dip in the Mbh 
is not as extreme. This is also true for the peaks, with the 
exception of the first peak of the later starburst/accretion 
events in SPHS at t ~ 1 Gyr. We have already mentioned 
how the second SPH starburst is more powerful than its 
SPHS equivalent, and this is seen also in the SFR, with the 
SPH peaking at a higher rate and falling off more slowly. 

From this point on our comparisons are made with the 
higher resolution set of runs. We outline a few specific differ- 
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Figure 1. Projected surface density plots for SPH-96-resO (left) and SPHS-96-resO (right), both with A^gas = 1.5 X lO'^, at i = 2 Gyr. 
There is a clear difference in the amount of structures present in the hot halo due to supernovae feedback - in SPH the shape of the 
overpressuriscd bubbles is more defined and a number of spherical clumps have condensed out of the ambient gas. In SPHS the gas 
distribution is more homogeneous, and no such clumps have formed. 



ences that we wish to highlight between the results obtained 
with the SPH and SPHS methods. 



3.1.3 Clumps vs. filaments 

Of particular interest is the behaviour of the two methods 
just after t = 1 Gyr. As we have already mentioned, the 
SPH run has started to see condensing clumps form out 
of the halo gas, which fall into the centre and form stars, 
giving rise to a powerful second starburst. As the bubble 
expands through the surrounding gas, there appears to be a 
sudden condensation of multiple dense clumps from the pre- 
viously low-density gas. These clumps scatter off each other 
as they fall to the centre, but remain seemingly protected 
from the ambient gas through which they travel. The clumps 
condense to the point that they reach the polytrope (refer 
to Sections 12.41 & 12. 5p and therefore undergo star forma- 
tion and subsequent feedback. The result is that the flow in 
the inner ~ 100 kpc becomes very disordered, with contin- 
ued asymmetric feedback events of varying strengths. The 
clumps that reach the central few kiloparsecs either add to 
or remain in orbit of the disc, although some are accreted. 

In the SPHS run, however, this condensation of cold, 
dense clumps does not occur. Rather, the structures that 
form are filamentary in nature, forming at the intersection 
between bubbles, with gas infalling through the filaments to 
grow the disc (or to be accreted by the SMBH). The very 
few spherical overdensities that do begin to form are quickly 
disrupted by travelling through the ambient gas before they 
can reach the densities required to lie on the polytrope. This 



quite striking comparison can be seen clearly in Figure O 
where the structure formation and mode of feeding is clearly 
different between the two methods. In Figures [3 &|8] (second 
from top, left-hand panels), we isolate the overdense regions 
through a series of cuts in the phase diagram. The differ- 
ences between the cold clumps (SPH) and filaments (SPHS) 
can clearly be seen. In Section 13.31 we discuss the reasons 
for these differences and demonstrate that the clumps are a 
numerical artifact. 



3.1.4 Disc formation 

Figure |3] shows the surface density projection of the gaseous 
disc that forms in the higher resolution SPH and SPHS runs. 
Both discs have a similar mass of ~ 2 X 10^ Mq , and are of 
a similar size, extending out to « 2.5 kpc. These similarities 
are understandable given the identical initial condition used 
for both runs. However, the formation of the disc is strongly 
influenced by the asymmetries present in the star formation 
and feedback events, which eject hot bubbles in a variety of 
directions into the gaseous halo. It is therefore not surprising 
that the discs in the SPH and SPHS runs are at different 
orientations, as the star formation and feedback histories 
of the two are somewhat different. The stochastic nature 
of the feeding to the central regions is so significant that 
neither disc is in the preferred x — y plane set by the angular 
momentum of the initial conditions. The two discs do in 
fact start out in the x — y plane as they are initially formed 
in the first starburst, but undergo precession and midplane 
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Figure 2. Projected surfaee density plots for SPH-96-resl (left) and SPHS-96-resl (right), both with A^gas = 7.5 X 10^, at t = 1.2 Gyr. 
While there is structure in both runs at this resolution, the mode of structure formation is very different, with the gas that has cooled 
to ~ 10* K condensing into a myriad of clumps in the SPH run but forming into two distinct filaments in the SPHS run. As a result, the 
manner in which the central gaseous disc grows is different; in SPH some of the clumps scatter off one another and fall into the centre 
to feed the disc, whereas in SPHS the disc grows as gas is funnelled down through the dense filaments that form due to the intersection 
of two or more supernovac-driven bubbles. 



rotation as they are torqued by subsequent infalling gas at 
a variety of different orientations. 

One important distinction between the discs in the two 
methods is the presence of the clumps orbiting outside the 
disc in the SPH case, which are not there in SPHS. This 
relates to the clumpy vs. filamentary mode of feeding the 
disc as discussed above in Section 13.1.31 There are clumps 
within the disc in the SPHS run, but they have formed out of 
the disc itself through gravitational instability, rather than 
being formed further out in the halo and migrating to the 
centre. 



3.1.5 Galaxy morphology and stellar content 

When we look in more detail at the kinematics of the gas 
and stars in the inner region, we find that the morphologies 
of both the stars and the gas actually differ by consider- 
ably more than a visual inspection of the disc (see Figure 
[3| would suggest. In particular, the SPH-96-resl run shows 
a larger amount of low angular momentum material within 
the central 8 kpc than is present in the SPHS equivalent. 
This can be seen clearly in Figure [l] which shows plots of 
the Jz/Jc ratio (often referred to as the 'disc/bulge' ratio), 
where Jz is the magnitude of the angular momentum vector 
along the rotation axis of the inner 2 kpc of each disc, and 
Jc is the magnitude of the angular momentum for a circular 
orbit at a radius r, namely Jc — rvc, where = GM{r)/r 
is the circular velocity from the enclosed mass of all species 



(gas + stars + dark matter -I- SMBH). This ratio was cal- 
culated for each star and for each gas particle and the his- 
tograms for both shown in the Figure (top left and right 
respectively). It is clear from these plots that the SPHS run 
shows a kinematically more pronounced disc than SPH, par- 
ticularly in the gas but also in the stars as well. The many 
cold clumps that are formed in SPH provide the inner 8 kpc 
with a greater amount of low angular momentum material 
(gaseous and stellar) as well as a greater randomness of or- 
bits as they fall in. The top right plot tells us further that 
the gas disc seen in SPH in Figure [3] is more disturbed and 
warped compared to the SPHS one, being strongly torqued 
out of its plane and heated by the infalling gas clumps at 
different orientations. We note that since it is this gas that 
tends to lie on the polytrope and is therefore 'star- forming', 
the stars that form here will likely inherit the Jz/Jc ratio 
seen for the gas in the top right plot. This means that to 
later times (e.g., 3 — 4 Gyr), the stellar kinematic disc/bulge 
ratio will likely become even more disparate between the 
SPH and SPHS simulations, with SPH favouring a bulge 
and SPHS favouring a disc. 

Also in Figure |4] we see a plot of the enclosed mass in 
stars with radius (bottom left) and enclosed mass of gas 
with radius (bottom right), showing that the stellar content 
of the galaxies in SPH and SPHS differ considerably outside 
of '-^ 2 kpc. The SPH run has a greater amount of mass in 
stars at all radii, and in particular further out in the halo, 
where the cold clumps are forming stars out to 20 — 30 kpc. 
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Figure 3. Projected surface density plots of tlie central disc that forms in the SPH-96-resl (left) and SPHS-96-resl (right) runs, at 
t = 2 Gyrs. Both discs are of a similar size, although the stochastic mode by which they have grown is very dependant on the random 
orientations of the supernovae feedback, and they are therefore at different orientations despite the same initial velocity field for the 
gaseous halo. The SPH disc is surrounded by multiple clumps that formed in the halo whereas the clumps present in the SPHS case have 
formed from the disc itself. 



In the gas the profiles are similar, but the SPHS run has a 
greater mass of gas at large radii in the hot halo. 

3.2 Features in p — T space 

In order to follow the relevant behaviour and understand 
evolution of the simulations better, we analyse a number 
of aspects of them in terms of their locations on a phase 
diagram p vs. T. We isolate a few key regions and follow the 
behaviour of the gas that inhabits these regions on the phase 
diagram. The corresponding plots for this section are Figures 
[7] & [SI in which the phase diagrams are plotted along with 
the projected surface density and projected temperature for 
the isolated phase regions. 

3.2.1 Isothermal contraction phase 

As the gas cools out of hydrostatic equilibrium, it 
experiences a variable cooling rate that is se t by 
the combined cooling curve of iKatz et al. fl996l ) and 
iMashchenko et al. (2008 ) used in our simulations. At T = 
lO'' K (the switchover between the two curves) the cooling 
becomes inefficient, and so gas tends to 'pile up' at around 
this temperature. We have isolated a region on the phase 
diagram (2nd from top) that corresponds to gas that has 
cooled to this temperature and is condensing, moving nearly 
isothermally across the marked region to higher densities 
and eventually onto the polytrope. This we refer to as an 
'isothermal contraction phase', and has particular relevance 



to the many spurious clumps that are seen in the SPH runs, 
as it marks the region on the phase diagram that they in- 
habit. Although still present in the higher resolution SPHS 
runs, this region is significantly less occupied, and does not 
extend as far to low densities as in the SPH case. 



3.2.2 Polytrope phase 

As mentioned in Sections 12.41 fc the diagonal line shown 
on the phase diagram corresponding to a polytropic EQS 
with an adiabatic index of 4/3 functions as a dynamic pres- 
sure floor, ensuring that the Jeans mass is always resolved in 
the gas. The polytrope 'fills in' from gas that either cools di- 
rectly onto it or from gas that reaches it through the isother- 
mal contraction phase discussed above. Once on the poly- 
trope, the gas can move along it to higher densities (and 
higher temperatures as per the EQS) and of course can move 
off it either by an increase in temperature for a given den- 
sity, or by a decrease in density for a given temperature. Gas 
can also leave the polytrope by being converted into stars, 
which is allowed to occur beyond the fixed threshold of 10^ 
atoms cm~^, denoted by the green dotted line. The presence 
of gas on the polytrope past this fixed threshold is therefore 
transitory - it fills in and subsequently disappears as stars 
form. 
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Figure 4. Plots of the kinematic 'disc/bulge' ratio at t = 2 Gyr for the inner 8 kpc in the stars (top left) and the gas (top right) together 
with the enclosed mass in stars (bottom left) and gas (bottom right) out to 100 kpc for the SPH-96-resl (black solid) and SPHS-96-resl 
(blue dashed) simulations. The SPHS run has a more disc-like morphology inside the central 8 kpc, with the Jz/Jc — 1 signal (a 'disc') 
being a factor of i?i 4 higher than SPH in the stars and a factor of si 5 higher than SPH in the gas. The SPH run shows a stronger peak 
at lower angular momentum with a Jz/Jc — signal (a 'bulge') that also corresponds to more randomised orbits than in SPHS. 



3.2.3 Hot bubble phase 

The filling in of the polytrope beyond the fixed density 
threshold for SF is often followed by an ejection event that 
pushes gas into the 'overpressurised bubble' section of the 
phase diagram. Within this region there we can identify both 
the p — T trend for a given bubble, as well as the evolu- 
tion of this trend with time, through analytic arguments. To 
start with, the similarity solution for a Sedov- Taylor blast 
wave due to an energy deposition S in a uniform medium 



(jSedov. I959I ) is given by; 

r{t) oc ( — ) (19) 

\ PISM / 

where pisM is the average (constant) density of the sur- 
rounding medium. We can define an overdensity parameter 
Ss = ps/pisM for the shocked gas, which under the assump- 
tion of a strong adiabatic shock, is a constant - given by 
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Figure 5. SFR and SMBH accretion histories for Afgas = 2 X 10^, with SPH (left) and SPHS (right). Sohd line is SFR, dotted is A/bh- 
The initial starbursts and accretion events start off identical, but in SPHS both the SFR and Af^h a^rc prolonged compared to SPH, with 
a secondary peak just as the main starburst is shutting off. The SPH run shows a second major starburst not present in the SPHS run. 
In both simulations the shape of the SMBH accretion rate curve largely follows that of the SFR, albeit with a small offset in time. 
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Figure 6. SFR and SMBH accretion histories for Afgas = 1 x lO'^, with SPH (left) and SPHS (right). Solid line is SFR, dotted is M^h- 
The initial starburst has a slightly higher SFR in SPHS but the initial accretion rate is higher in SPH. Both rates then drop considerably 
for the next ~ 1 Gyr with the SFR dropping to almost zero and the accretion rate by several orders of magnitude. During this period 
the SPH accretion rate remains higher than in SPHS by a factor of ~ 10. The second major starburst occurs later in SPHS than in SPH, 
has a lower SFR and drops off faster. However, in SPHS there are two subsequent major starbursts rather than just one in SPH. Again, 
the SMBH accretion rate approximately follows the shape of the SFR, with a small offset. 



Ss — 4. The shock velocity, Vs = dR/dt is therefore: 

o or 
which we can write as: 



(20) 



The post-shock temperature for an adiabatic shock with ve- 
locity Va is: 



which gives us: 



Ta cc 



(22) 



rpl/2 -3/2 -1/2 
Vs (X E ' r Pa 



(21) 



Ta OC Er-^'pJ^ 



(23) 
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Figure 7. Projected surface density plots (left), phase diagrams (middle) and projected temperature plots (right) in SPH-96-resl at 
t = 1.2 Gyr. The white dot at (0,0,0) marks the SMBH. The blue dotted line indicates the cooling floor, while the green and brown 
dotted lines mark the fixed SF threshold and polytropic pressure floor respectively. The two thicker diagonal lines indicate the standard 
Sedov- Taylor solution (magenta) and the post Sedov- Taylor adiabatic expansion (black) that fit the orientation of the feedback peaks and 
their evolution in time respectively. The red outlined boxes indicate the regions of the phase diagram being plotted. Contours are binned 
logarithmically, starting at the minimum value for p or T and increasing with a factor of dp/p = dT/T = 0.1. The colours correspond to 
the number of gas particles, with cyan representing individual particles and subsequent colours going up in factors of 5. In this simulation 
we see a great deal of structure, in the form of multiple overdense clumps, distributed over a large part of the computational domain and 
lying within the 'isothermal contraction' region of the phase diagram. We also see a broadened spike along the Sedov- Taylor gradient 
that corresDonds to manv different suoernova feedback events. 
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Figure 8. Projected surface density plots (left), phase diagrams (middle) and projected temperature plots (right) in SPHS-96-resl at 
t = 1.2 Gyr. Linestyles etc. arc as per Figure [7] There are striking differences in the structures that have formed compared to the SPH 
case (Figure O - in the latter we saw multiple clumps whereas here we see one or two filaments. Again, the structure occupies the 
isothermal contraction region of the phase diagram, with the gas flowing down the filaments to feed the disc and form stars. Similar 
to the SPH case, wc sec a broadened spike along the Sedov- Taylor gradient corresponding to multiple feedback events, but reduced in 
magnitude. 
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The trend T (x p^^ is shown on the top phase plot (dotted 
magenta hne) in Figures [7] & [8] The shape of this feature 
remains constant but evolves along the phase diagram as 
the bubbles expand. Such evolution can be described ana- 
lytically by a post Sedov- Taylor solution for adiabatically- 
expanding hot gas. If we assume that the expansion is suf- 
ficiently fast that radiative cooling can be neglected, the 
equation of state is P oc p^, with 7 = 5/3. We can write this 
in terms of the sound speed by noting that: 



and therefore that: 

(25) 

which we equate to obtain: 

T (X p^'^ (26) 

for the time evolution of the bubble. This trend is also plot- 
ted in the top phase diagram (dotted black line) in Figure 
[7]&|8] 

Even though the analysis described above corresponds 
to the solution for a uniform medium, it fits the supernovae 
ejecta in our simulations remarkably weljf]. The reason for 
this is that the two departures from the standard Sedov- 
Taylor solution that are present in our model exert compet- 
ing effects on the p — T profile of our hot bubbles. Firstly, we 
have a non-uniform medium, where p falls off according to 
equation [TT] This has the effect of reducing the velocity fall- 
off as the shell propagates and sweeps up mass. Secondly, 
we have an inwardly-directed velocity field as the gas cools 
and inflows under the influence of the gravitational poten- 
tial, which has the effect of enhancing the velocity fall-off. 
The infall velocity arises self-consistently from the potential, 
which is a function of the density profile p{r). The two de- 
partures from standard Sedov- Taylor therefore balance out, 
and the majority of the supernovae-driven bubbles behave 
according to the analytical arguments above (with small ad- 
ditional departures due to radiative cooling and the modi- 
fication of the density profile due to previous supernovae 
explosions) . 

3.3 Numerical clump formation 

So far, it is clear that the presence of the overdense clumps is 
a feature of the SPH method, but the question remains as to 
what causes them. The fact that an identical simulation run 
with SPHS (CS kernel with 96 neighbours) does not yield 
these structures suggests that the prevention of multivalued 
fluid quantities plays a role in avoiding their formation, as 
this is the main difference between the methods when the 
kernel and neighbour number a re not sufficient to have im- 
proved force accuracy (refer to iRead fc Havfield. 20121 ) . In 
this section, we show that this is indeed the case; it is the 
removal of pressure blips in an otherwise smooth flow that 
prevents the condensation of the clumps. 

To iden tify bound gaseous clumps, we ran the Am iga 
halo finder (|Gill et al.. 2004KnoUmann fc Knebe. 20091 ) on 

^ to see the time evolution clearly the reader is directed to 
|http : //www . phys ■ ethz . ch/ - ahobbs/movles ■ html | 



each simulation output. We focus here on the most gas- 
rich clump in the SPH run. Tagging the clump particles, we 
traced the evolution of the clump back in time to its initial 
formation. This particular clump forms from the merging of 
three distinct regions of geis, one of which has a significantly 
higher density and lower temperature than the other two. 
As the three regions merge, a small density peak forms with 
a corresponding entropy dip caused by gas cooling in the 
peak. The flow at this point is strongly shearing and the 
density peak rapidly shears away. By contrast, however, the 
entropy dip remains due to the lack of multiphase mixing in 
SPH. The presence of an entropy dip with no corresponding 
density peak drives a pressure dip at the centre of the merg- 
ing gas. This can be seen on the left-hand side of Figure 
[51 where the entropy and pressure dips are marked by the 
red and magenta circles, respectively. The central pressure 
dip drives an inward collapse and the formation of a bound 
clump (see Figure [TCT)) . 

In order to be sure that this multivalued pressure prob- 
lem is not present in the new SPHS method, we took the 
starting conditions from the SPH run at t = 0.4 Gyr - just 
before the 'clump gas' began to merge - and ran it with 
SPHS (CS kernel, 96 neighbours). Once again, we tracked 
the evolution of the gas which, in the SPH case, ended up 
in the clump. The initial merging occurred in exactly the 
same way, but as the density peak sheared away, so too did 
the dip in the entropy. As a result, there was no pressure 
dip and therefore no initial seed for collapse. This can be 
seen on the right-hand side of Figure |9l where the features 
marked by the circles on the left-hand side are not present. 

Figure [To] shows the state of the same 'clump-identified' 
gas in the SPH and the SPHS cases. The former, as a re- 
sult of the 'pressure-dip seed', has collapsed to form a near- 
spherical clump with multiple overdensities. In the SPHS 
run, however, the gas has been drawn out into a filament, or 
stream, having not been artificially pulled inwards due to a 
central dip in the gas pressure. 

The filament/stream seen in Figure [TOl for the SPHS 
evolution of the merging clump shows individual density 
peaks forming along the central axis - however, at this res- 
olution these are dissimilar structures from the clumps seen 
in SPH, as they are transient, being mixed in with the sur- 
rounding gas before they can form stars. In the next Section 
we discuss the formation of bound clumps within the fila- 
ment at higher resolution, which form due to the density 
enhancement created as two or more evacuated bubbles in- 
tersect. 



4 'RESOLVED' CLUMP FORMATION WITH 
SPHS: FRAGMENTATION OF OVERDENSE 
FILAMENTS 

We now move on from the comparison between the two 
methods at identical spatial resolution to a simulation that 
employs the full SPHS algorithm: in addition to the dissi- 
pation for advected fluid quantities, we also have improved 
force accuracy through the use of the HOCT4 kernel with 
442 neighbours. We perform this simulation at significantly 
higher resolution than our fiducial comparison runs, reach- 
ing a gas particle mass of a few xlO* Mq (see Table 1). 

The evolution of the higher resolution simulation is 
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Figure 9. Plots of density (top), pressure (middle) and entropy (bottom) for the most gas-rich overdense clump in the SPH-96-resl run, 
at t = 0.63 Gyr. The gas that formed the clump was identified and tracked back earlier in the simulation to before it formed. On the 
left-hand side are the plots showing the evolution of this particular clump with the SPH method, while on the right-hand side the plots 
show the evolution of the clump with the SPHS method (although starting from the same SPH-96-resl snapshot from which the clump 
was taken — at t = 0.4 Gyr). Each property (p, P, A) is plotted in the centre of mass frame of the clump gas at a time in its evolution 
where a density peak that occurred just previously has been smoothed out. This density peak had a corresponding entropy dip at the 
same location. When evolved with SPHS, the smoothing of the density peak coincides with a diffusion of the entropy, removing the dip 
and allowing the pressures to remain smooth; however, when exactly the same initial condition is evolved with SPH the entropy dip 
remains (red circle), driving an equivalent dip in the pressures (magenta circle). This central pressure dip drives the contraction of the 
gas and the subsequent formation of the clump. For the full evolution and a clear picture of how this occurs the reader is directed to the 
corresponding movie at |http: //www.phys ■ ethz. ch/- ahobbs/movles .htmll 
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Figure 10. Visualisation of the numerical experiment discussed in Section [3.3l and Figure[9] The initial condition (top plot) at f = 0.4 Gyr 
is identical (taken from the SPH simulation of an earlier state of the gas in one of the bound clumps) but was evolved both with SPH 
(bottom left) and with SPHS (bottom right) to a time of f = 1 Gyr. The differences are striking; the SPH evolution of the gas has 
formed a near spherical, dense clump while the SPHS evolution shows a stretched out filament of lower density. The point at which the 
two diverge is shown in Figure |9l where a pressure discontinuity at the centre of the merging gas in the SPH case causes an artificial 
contraction that leads to the formation of the clump; this feature is not present in the SPHS case. 
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globally similar to those at lower resolution - an initial cen- 
tral starburst sends a strong Sedov- Taylor bubble into the 
surrounding gas, which shuts off star formation and reduces 
SMBH accretion until such time as the gas makes its way 
back to the central regions. The subsequent starbursts are 
then highly asymmetric, with multiple cavities created from 
individual supernovae feedback events. Again, there is no 
sign of the many artificial clumps condensing from the am- 
bient halo gas, and the structure that forms as the gas cools 
is filamentary in nature, caused by the merging of the outer 
walls of two or more of the feedback-driven bubbles. 

Here, however, we are in a new regime whereby bound 
clumps are in fact able to form from the hot halo, but 
only within the filament (s). The latter structure provides a 
strong density enhancement within which collapse and frag- 
men tation can occur through non-linear the rmal instabil- 
ity (jjoung et al.. 201 ^Fernandez et al.. 2012l l. rather than 
through a numerical pressure discontinuity as was the case 
for the clumps in SPH. Figure[TT]shows the state of the inner 
30kpc at a time when the filaments have recently formed. 
We used AHF to identify the locations of any self-bound 
clumps; these are shown in blue. One or two have already 
collapsed and are forming stars. As can clearly be seen from 
the Figure, the only locations where clumps arise are in the 
streams and in the disc at the centre. The overdensity here is 
~ 10 — 100, which is consistent with the instability t hreshold 
for a non- linear perturbation in ljoung et al. 

The gas flowing through the filaments ends up con- 
tributing to the disc, as do the clumps that form from the 
fragmentation of the filament. Some of these, however, form 
stars before they can reach the disc, and therefore end up 
as orbiting stellar clusters, with masses of ~ IO'^Mq. The 
stellar clusters, while they are still star-forming, contribute 
to the continued feedback events and can therefore aid the 
formation of subsequent streams that may feed the disc. 

In order to quantify the ability of the cold gas in the 
filaments to grow the disc, in Figure [12] we plot the gas 
mass contained in particular phase regions (the same regions 
marked in Figures [7| & [8] but for our highest resolution full 
SPHS run). We also plot the rate of increase or decrease in 
mass within this region to measure the feeding rate through 
the filaments. The main period of interest is after the sec- 
ond starburst, and we see from Figure [12] (top plot) that the 
streams (and the clumps formed in the streams) are feeding 
the central disc at a rate of ~ 1 Mq yr~^ for the first Gyr af- 
ter the starburst. This rate is gradually declining over time, 
and by t = 2Gyrs has dropped to 0.1 yr^^. Looking 
at the middle plot, we see further that the gas on the poly- 
trope (most of which is in the disc) exhibits a similar rate 
of decrease in mass; this can be explained by referring to 
the SFR (bottom) plot in Figure [TTI which shows that the 
feeding rate of the disc by the filaments is approximately 
matched by the star formation rate, since most of the star 
formation is occurring in the disc at this time. 

The 3rd (bottom) plot of Figure [12] shows the amount 
(and rate of increase/decrease) of gas in the region of the 
phase diagram that we have identified as corresponding to 
the gas recently ejected from supernovae explosions. We see 
here that this is by far the dominant repository for the gas 
that has already fallen into the central regions - the mass 
contained in stars, disc, filaments and clumps combined adds 
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Figure 11. Projected surface density plot of the SPHS-442-res2 
run at f = 0.73 Gyr. As in the lower resolution case, we see the 
formation of one or more cold filaments that have formed out of 
the surrounding e result of the interaction of supernovae- 

driven hot bubbles. The filaments are funnelling gas down onto 
the disc at the centre of the computational domain. The progen- 
itors of gas-rich clumps that are formed in the filaments and in 
the disc are marked in blue. 



up to ^ 1/lOth of the mass that has been ejected into the 
hot halo, the latter being on average ~ 1-2 x 10^" Mq . 



5 DISCUSSION 

We find a clear difference in how cold gas condenses from 
a hot halo in SPHS versus 'classic' SPH. The formation of 
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~ hundreds of cold clumps from relatively homogeneous re- 
gions of the hot halo in SPH owes to a numerical thermal 
instability. As hot supernovae-driven bubbles collide in the 
halo, the gas between the bubbles is compressed to high 
density and cools. The flow is highly shearing and the over- 
density rapidly shears away. Due to the lack of multiphase 
mixing, however, the gas remains cold leading to a pressure 
dip that seeds the formation of a dense clump of gas. By 
contrast, in SPHS the shearing gas mixes both in density 
and entropy and the clumps do not form; instead the gas 
forms into cold filaments that feed the disc. In our high- 
est resolution SPHS run, these filaments break up to form 
bound clumps. This fragmentation is physical and owes to a 
non-linea r instability caused by the ^ 10 overdensity in the 
filament (|Joung et al.. 20121 ). 

iKaufmann et al. invoke numerical noise as a po- 

tential seed for small-scale fluctuations that could lead to 
runaway collapse. Indeed, we find that at least for some of 
the clumps that we see forming in SPH, the seeds are per- 
haps one or two particles with significantly lower entropies 
(and initially higher densities) than the gas they end up in as 
regions of the simulation merge. These lower entropies can 
lead to pressure dips (P oc Ap"* , where A = A{s), a, function 
of the specific entropy s) when the densities are smoothed 
across the neighbours, as in SPH the entropy of a particle is 
'protected' and is not shared. In SPHS this 'entropy noise' 
is effectively removed via the sub-grid diffusion. 

There may be a situation where the formation of such 
clumps from thermal instability is physical, as indeed is 
suggested bv IKaufmann et al. (20091 ) and explored analyt- 
ically by iBinnev et al. f2009l ). which is when the entropy 
profile across the region is nearly completely flat and the 
cooling timescale can become smaller than the oscillation 
timescale for small, isotropic perturbations. However, such 
an entropy proflle is highly unlikely for a galaxy formation 
simulation (both cosmological and non-cosmological) and 
indeed is certainly not the case in our simulations, where 
the entropy gradient is quite steep both initially and af- 
ter mul tiple supernovae fee dback events. Even in this sit- 
uation, iBinnev et al. (20091 ) find that thermal conduction 
suppresses linear thermal instabilities through the damping 
of small-scale p erturbations, e ven for very small fractions of 
Spitzer's value (|Spitzer. 19621 ). 

Our finding that at sufficient resolution the SPHS 
simulations show resolved clumps forming in the over- 
dense fllament(s) that feed the disc has particular r ele- 
vance to the recent w ork bv iFraternali fc Binnev (20081 ) & 
iMarinacci et al. (201ll ). where the authors suggest a model 
of gas entrainment from the hot halo through the motion 
of feedback-driven galactic fountain clouds. These colder 
clouds mix with the coronal gas and invoke a significant 
amount of mass transfer from the hot phase to the cold 
phase, before sinking back again to the galaxy. The key 
concept here is that the cooling of gas from the halo is 
dependant on the presence of star-forming gas deeper in- 
side the potential well. Our picture is slightly different, 
but still within the same vein, as it is the interaction of 
the walls of multiple supernovae-driven bubbles that cause 
the condensing of the filaments, which in turn may break 
up and form cold clum ps through non-linear instability as 
per ljoung et al." (20121 ). Naturally, therefore, we require that 



star formation occurs in the disc in order for cold gas to be 
able to accrete - a for m of positive feedback. 

Observations by iBoomsma et al. (2005h lend support 
to this picture, where X-ray emission from hot gas is 
seen above the plane of the nearby spiral galaxy NGC253, 
along with an overdense filament-like structure detected in 
HI at the edge of the hot region, that stretches out for 
approximately 12 kpc. These features are associated with 
a significant amount of star formation activity and in- 
deed a recent starbu rst and superwind from the galaxy 
(|Heckman et al. 199(]| ). 

In our highest resolution fuU-SPHS simulation, we 
find that the fragmenting filaments formed from over- 
lapping supernovae-driven bubbles feed the disc (and 
subsequently, star formation in the disc) at a rate of 
~ IM0 yr~^, similar to what is required by ob- 
servations to fuel star formation in th e Milky Way 
(|Noh fc Scalo. 199ClRocha-Pinto et al.. 200Cl ). We also find 
that the gas mass ejected from supernovae explosions re- 
mains near-constant around 1-2 x 10^" M0 for the duration 
of the simulation (2 Gyrs). If one assumes that the inner hot 
halo is composed primarily of gas from supernovae feedback 
processes then this is in excellent agreement with the mass 
of th e hot halo as determined by pulsar dispersion mea sures 
(e.g.. I Anderson fc Bregman. 201(IGaensler et al.. 20081 ). We 
emphasise, however, that our goal in this paper is to study 
the physics of thermal instabilities in cooling haloes rather 
than to build a faithful model of Milky Way-mass disc galax- 
ies. Our results are likely sensitive to our sub-grid model 
parameters, to our resolution, and to our initial conditions 
that are idealised and non-cosmological. Thus, while we can 
identify an important thermal instability in the halo that 
can in principle supply significant cold gas to the disc, it is 
unclear how important it is with respect to other mecha- 
nisms like gas rich mergers, low-redshift cold streams, and 
gas recycling from stars in the disc. We will explore this in 
future work. 

Finally, a further aspect of the comparison between 
the SPH and SPHS runs at identical resolution that is 
particularly interesting is the difference in the kinematic 
morphologies of the galaxies due to the different cooling 
mechanisms, i.e., clumps vs. filaments. In the SPH case the 
cold clumps contributed to a stronger galactic bulge-like 
feature with a greater amount of material at low angular 
momentum and on random orbits, while in SPHS the 
growth of the disc from the cold streams led to a stronger 
disc feature that was less warped and undergoing a lesser 
torque from infalling material at different orientations. 
This result naturally has implications for the so-called 
'angular momentum problem' in galaxy formation, namely 
the presence of bulges that are over-dominant and discs 

that are under-dominant compared to observations (see, e.g., 

iMaller fc Dekel, 200 jBurkert fc D'Onghia. 2004Governato et al.. 2010l ). 



6 CONCLUSIONS 

We have presented simulations of a cooling gaseous halo in 
a Milky Way mass galaxy, using this particular problem to 
perform the first scientific investigation with a new hydrody- 
namics code, SPHS. We have compared the results obtained 
(at identical spatial resolution) with that of the standard 
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('classic') SPH method employed in the literature, finding 
significant differences in the mode of gas cooling in the halo 
and subsequently the mode of disc feeding. The puzzle of 
the many cold clumps seen forming from the halo in many 
SPH simulations of galaxy formation is attributed to a nu- 
merical inability driven by unresolved mixing of different 
gas phases. We demonstrate both with our full simulations 
(Section I3.1|l and -with a more idealised test (Section I3.3|l 
that the removal of pressure blips in an otherwise smooth 
flow prevents the formation of the clumps, leading instead 
to the formation of cold filaments that feed the disc. The 
resulting galaxies in the SPH and SPHS simulations differ 
greatly in their morphology, gas phase diagrams, and stellar 
and gaseous disc/bulge ratio. 

We have explored in more detail the mode of disc feed- 
ing seen in our SPHS simulations, going to higher resolution 
and employing a kernel that gives improved force accuracy. 
We find a new way of bringing cold gas to the galactic disc; 
namely, the fragmentation and collapse through non-linear 
thermal instability of filament(s) formed at the intersection 
of supernovae-driven bubbles. The feeding rate of cold gas 
(T ^ 10'' K) to the disc is found to be approximately a solar 
mass per year, which suggests this is a promising model for 
fuelling late-time star formation in real spiral galaxies. 

We emphasise that our focus in this paper was on un- 
derstanding what drives thermal instabilities in hot halo gas, 
and in particular performing a comparison between the 'clas- 
sic' SPH and the SPHS numerical methods. Our numerical 
experiments are idealised and do not present a complete 
picture of galaxy formation. Nonetheless, by employing a 
hydrodynamics method that resolves the mixing of different 
gas phases, we find a novel mode of cold gas accretion and 
disc growth that may be very relevant for galaxy formation. 
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Figure 12. Plots of the gas mass contained within a particular 
region of the phase diagram (black) along with the rate of change 
of the gas mass into (blue; an 'inflow') and out of (red; an 'out- 
flow') the phase region. The regions are as per Figure[8]but in this 
case for the SPHS-442-res2 run. The 'inflow' and 'outflow' rates 
arc plotted in bins of 80 Myr. The plots correspond to; (i) the 
gas in the 'isothermal contraction' region (top) - this comprises 
the temporary disc created in the initial starburst and then later 
the filaments that form and grow a new disc; (ii) the gas in the 
'polytropc' region (middle) - this comprises the star forming gas 
i.e., the clumps that form in the disc and the filaments; and (iii) 
the gas in the 'hot bubble' region, corresponding to gas that has 
been ejected into the hot halo by supernova feedback. 
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